#!/bin/bash
set -e

while getopts  ":i:p:" opts
do
        case  $opts  in
		i)
            interval=$OPTARG       
                ;;
		p)
			out_prefix=$OPTARG
                ;;
		\?)
			echo `basename $0` [-p out_prefix] [-i interval.bed] 
			exit 1 
				;;
        esac
done
shift $(($OPTIND - 1))

# if [ -z $1 ]; then
	# echo `basename $0` [-p out_prefix] [-i interval.bed] 
	# exit 1
# fi



. /mnt/ilustre/app/medical/tools/.var


# genome_name=b37.fa
echo using $genome_name as the reference genome
# genome_assembly=b37
echo genome assembly is: $genome_assembly
# dbsnp_version=138
echo dbsnp version is: $dbsnp_version

# data_thread_num=8
# cpu_thread_num=4

# java_memory=16g
echo java memory: $java_memory

# snpeff_db_version=GRCh37.75
echo snpeff database: $snpeff_db_version


vcf_path=${data_path}/vcf/gatk/

echo ref genome is: $ref_genome


reads_seq_1=$1
reads_seq_2=$2

echo reads are: $reads_seq_1 and $reads_seq_2


if [ -z "$interval" ]; then

echo
echo
echo gatk BaseRecalibrator
java -Xmx$java_memory -jar $gatk \
	-T BaseRecalibrator \
	-R $ref_genome \
	-I $out_prefix.realn.bam \
	-nct $cpu_thread_num \
	-knownSites $data_path/ncbi/dbsnp/All_20150605.vcf \
	-knownSites ${vcf_path}Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf \
	-knownSites ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-o $out_prefix.realn.recal_1.table

else

echo
echo
echo gatk BaseRecalibrator
java -Xmx$java_memory -jar $gatk \
	-T BaseRecalibrator \
	-R $ref_genome \
	-I $out_prefix.realn.bam \
	-L $interval \
	-nct $cpu_thread_num \
	-knownSites $data_path/ncbi/dbsnp/All_20150605.vcf \
	-knownSites ${vcf_path}Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf \
	-knownSites ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-o $out_prefix.realn.recal_1.table
	# -rf BadCigar
	
fi

echo
echo
echo gatk PrintReads
java -Xmx$java_memory -jar $gatk \
	-T PrintReads \
	-R $ref_genome \
	-I $out_prefix.realn.bam \
	-nct $cpu_thread_num \
	-BQSR $out_prefix.realn.recal_1.table \
	-o $out_prefix.realn.recal.bam

	
. $cmd_done
# echo
# echo
# echo mark primer
# ${script_path}/mark_sam_qual0_start_end.sh aln.$out_prefix.realn.recal.bam mark.bam

# mv $out_prefix.realn.recal.bam $out_prefix.realn.recal0.bam
# ln -s $out_prefix.realn.recal.bai $out_prefix.realn.recal0.bai
# mv mark.bam $out_prefix.realn.recal.bam
